Pairing and Vortex Lattices for Interacting Fermions in Optical Lattices with Large 

Magnetic Field 



Hui Zhai, 1 Q R- 0. Umucalilar, 2 and M. 6. Oktel 2 

1 Institute for Advanced Study, Tsinghua University, Beijing, 100084, China 
2 Department of Physics, Bilkent University, 06800, Ankara, Turkey 
(Dated: April 14, 2010) 

We study the structure of pairing order parameter for spin- 1/2 fermions with attractive inter- 
actions in a square lattice under a uniform magnetic field. Because the magnetic translation sym- 
metry gives a unique degeneracy in the single-particle spectrum, the pair wave function has both 
zero and finite momentum components co-existing, and their relative phases are determined by 
a self- consistent mean-field theory. We present a microscopic calculation that can determine the 
vortex lattice structure in the superfluid phase for different flux densities. Phase transition from a 
Hofstadter insulator to a superfluid phase is also discussed. 



Optical lattices and synthetic magnetic fields are two 
of major tools to create strongly interacting many-body 
systems in cold atoms [HQ. In conventional solid state 
materials, accessible magnetic flux per unit cell ub is 
very small, ne C 1, even for the strongest magnetic field 
attainable in laboratory (< 45T). Hence, as in most con- 
ventional metals, the electron density n is several orders 
larger than ub that the magnetic field can be treated 
semi-classically; or as in the two-dimensional electron 
gases, n ~ wb < 1, the density is so low that only the 
bottom of an electron band is populated, and the effec- 
tive mass approximation is sufficient to account for the 
lattice effect. In cold atom systems, because the mag- 
netic field is synthetically generated by rotation [3| or by 
engineering atom-light interactions 0, [H[ , and the lattice 
spacing is of the order of half a micron, one can access 
the regime n ~ ub ~ 1, where both the lattice and the 
magnetic field should be treated on an equal footing and 
in a quantum- mechanical manner. Consequently, such a 
system exhibits the famous Hofstadter butterfly single- 
particle spectrum [6|. 

For neutral atoms in lattices, the interaction is dom- 
inated by on-site interactions as in the Hubbard model. 
Hereafter, we shall refer to the model describing inter- 
acting cold atoms in optical lattices with large magnetic 
field as the Hofstadter-Hubbard (HH) model. Recently, 
many works have focused on the bosonic HH model [7j, 
which reveal a number of interesting phenomena, includ- 
ing vortex lattice states and possible fractional quantum 
Hall states. However, so far little attention has been paid 
to the fermionic HH model. 

The subject of this letter is the properties of the paired 
superfluid phase in the fermionic HH model with attrac- 
tive interactions. For tlb ~ 1, the pairing problem differs 
from type-II superconductors in a fundamental way. In 
type-II superconductors the separation between the vor- 
tices is much larger than the size of Cooper pairs, hence 
one can locally apply the BCS scenario to define a local 
order parameter A(r), and understand the vortex lat- 
tice by coupling this "coarse grained" order parameter 



to the magnetic field. In the HH model considered here, 
magnetic field modifies the single-particle dispersion in 
an important way. Despite a strong magnetic field, there 
is always a well-defined Fermi surface and Bloch states 
in the magnetic Brillouin zone (MBZ) in the Hofstadter 
model. Therefore, with attractive interaction BCS pair- 
ing always occurs as an instability of the Fermi liquid. 
This enables us to reach the regime where the pair size 
is comparable to the distance between vortices, hence, 
any discussion of pairing must include the effect of the 
magnetic field at the microscopic level. We shall show 
that such a microscopic theory requires the definition of 
an order parameter with multiple components, and will 
discuss how this order parameter naturally describes the 
configuration of vortices. The main points of our analysis 
are highlighted as follows. 

(1) : We first review that for ub = p/q, where p and 
q are co-prime integers, each single particle state in the 
Hofstadter spectrum is o-fold degenerate due to magnetic 
translation symmetry [H, 0|. This degeneracy enforces 
that a comprehensive formulation of BCS theory in this 
case must contain Cooper pairs with both zero and a set 
of finite momenta, and treat them on an equal footing. 

(2) : We show that the magnetic translation symmetry 
also imposes relations between pairing order parameters 
of different momentum. These relations are verified nu- 
merically by self-consistently solving the BCS mean- field 
Hamiltonian. 

(3) : The relative phases between different pairing or- 
der parameters determined from self-consistent solutions 
can also be understood from a more intuitive and simpler 
Ginzburg-Landau argument. 

(4) : We determine the structure of vortices in the su- 
perfluid ground state using the information from (2). 
The unit cell of the superfluid phase is enlarged to q x q, 
whose symmetry is lower than that of the original Hamil- 
tonian. Hence, the superfluid ground state has discrete 
degeneracy, related to the symmetry of the vortex lattice. 

(5) : For certain fermion densities, a critical interaction 
strength is predicted for a quantum phase transition from 
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a Hofstadter insulator to a superfluid phase. 

The Model: We consider a two-component Fermi gas 
in a two-dimensional optical lattice potential so that an 
s— band tight binding model accurately describes the dy- 
namics. Both components are coupled to the same gauge 
field A = (0,px/q) in the Landau gauge. Note that there 
is no Zeeman shift associated with a synthetic magnetic 
field. The single-particle Hamiltonian is given by 

H = -tJ2 (e iW "'4<*<r + h.c.) , (1) 

where i = (i x ,i y ) labels the lattice sites, and (W) rep- 
resents all the nearest neighboring bonds. In the Lan- 
dau gauge, Aa> = pi x /q if i — i' is along ^/-direction and 
Aa> = if i — i' is along x-direction. Let T x and T y be 
magnetic translations of one lattice spacing along x— and 
y— direction (see Ref. [!, [sj for definition). The Hamil- 
tonian Hq commutes with both T x and T yi but these two 
operators do not commute T x T y — ex.p{i2irp / q}T y T x . 
One can choose the set of commuting operators as i7o,Ty, 
and (T x ) q = T qxi in effect, enlarging the unit cell in real 
space to contain q sites in x-direction. Thus, MBZ be- 
comes k x C [— 7r/q,7r/q) and k y C [— 7r,7r). Denoting the 
common eigenstate by i) n ,k x ,k y and the eigen-energy by 
tnk x k y , the magnetic Bloch theorem yields T q£ ip n ,k x ,k y = 
exp{iqk x }il) nykx}ky , T y ijj n , kx , ky = ex.p{ik y }i/j n , kx , ky , where 
n is the band index. Thus, for Tb, I = — 1, 

Ti x i/Jn,k x ,k y is a degenerate eigenstate of ^ n ,k x ,k y , and 
since T y (Ti x i/j n ^ x ^ ky ) exp{i(ky + 2nlp/q)}(Ti£^ n , kx ,k y ), 
we have the following properties [8|, |9[ 

1pn,k x ,k y (x + 1\V) OC ^n,k x ,k y +27rlp/q(x,y), 

£n,k x ,k y ~ £n,k x ,k y +27rlp/q' (2) 

For p/q = 1/3, the spectrum and Fermi surface shown in 
Fig. nja,b) clearly display a three-fold degeneracy. 

In addition to Ho, we consider the on-site interaction 
between different spin components 

The Hamiltonian for the HH model discussed below is 
then given by # H H = H + H int . 

Generalized BCS Theory: We start by diagonalizing 
the no n- interacting Hamiltonian Ho. First we relabel the 
sites to reflect the enlargement of the unit cell. For a site 
i = {i x ,i y ), we let i x — j x q + /?, where j x is an integer la- 
belling the magnetic unit cell and /3 = 0, . . . , q — 1 denotes 
the q- inequivalent sites within each cell. So, magnetic 
unit cells are uniquely labelled by j x and j y = i y . With 
this notation we identify Cjp a = Ci a . Fourier transfor- 
mation in the variable j yields Ck/3<j> where k is limited 
inside the MBZ. We now define a new set of operators 
4na through c k/5a = £) n (k)dkna, under which H be- 
comes diagonalized as 

H = ^ 6 ™ k 4n<An<x, (4) 
nkcr 
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FIG. 1: (Color online) (a) Three magnetic bands for p/q = 
1/3. (b) The Fermi surface of a half- filled (black solid line) 
and slightly away from half- filled (red dashed line) system, (c) 
and (d): the band dispersion along /^-direction with k y — 
(c) and along k y direction with k x = (d). Various possible 
pairings included in our BCS theory are also illustrated in (c) 
and (d), which include intra- and inter-band pairing (c), and 
pairing with non-zero center-of-mass momentum (d). 

Note that diagonalization of Ho is equivalent to solving 
Harper's equation [6], and ^2(k) is the /3th component of 
the nth eigenvector of Harper's equation at wave vector k 

• ^kn * s ^ ne operator that creates a particle in the nth 
magnetic sub-band at wavevector k. As an example, e n k 
is plotted in Fig. [TJa) for p/q = 1/3. In terms of dk na , 
Hi n t becomes 

#int = t/EEE E 

j3 Q k,k' ni ,...,n 4 

<?n§ + k ^n§ - wf - k x 4 (§ + k ') >< 

4 +k)ni)t 4_ k , n2 /§-k' )n3 ,; d #+k',n 4) t' ( 5 ) 

where the momentum sum is restricted to the MBZ. We 
should focus on the "on-shell" Cooper processes with 
k' = — k. Importantly, due to the (/-fold degeneracy, we 
not only consider Q = terms in Eq. (j5]), but also need 
to consider all the terms with Q = (0, 2-Klp/q), where 

1 = 0, . . . , q— 1, since — k and k + Q also have the same ki- 
netic energy. Consequently, non-zero center-of-mass mo- 
mentum pairing needs to be included as well. Besides, 
intra-band Cooper pairs have a non-vanishing coupling 
to the inter-band Cooper pairs. For instance, in Eq. (j5]), 
if rt\ = ri2 but ns ^ n±, the interaction coefficient is non- 
zero. Hence, intra-band pairing must induce inter-band 
pairing. All these pairing scenario under consideration 
are schematically illustrated in Fig. DJc-d), and a com- 
prehensive BCS theory in this problem must treat all 
these possibilities on an equal footing. Therefore, we in- 
troduce totally q 2 order parameters A^ = (Aq, . . . , A^_ x ) 
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^iven by 



n,n',k ^ ^ 



(d. 



(6) 



where = 0, . . . , g — 1. The site index /3 denotes 
g inequivalent sites along x-direction of each magnetic 
unit cell, and the index I represents the center-of-mass 
momentum of the pair Q = (0, 27rZp/g), which repre- 
sents the order parameter modulation along ^-direction. 
For instance, for p/q = 1/3, there are three differ- 
ent center-of-mass momenta, which are Qi=o = (0,0), 
Q i=1 = (0,2tt/3) and Q i=2 = (0,4tt/3). With AJ3, the 
mean-field Hamiltonian becomes 

nka 1,(3 I n,n' ,k 

(A^(k+ (-k+ §)< +§in /_ k+?i „, 4 



+h.c.) 



(7) 



The real space order parameter for site i = (i x ^y) is 
given by 



^(mod q) t 



i27Tlpi y /q 



(8) 



therefore the unit cell in the superfluid phase is enlarged 
to q x q in real space (see Fig. [2]). 

Solution to BCS Theory: We start with q 2 random 
complex numbers as initial and iteratively solve 
the BCS mean- field Hamiltonian [Eq. (|IJ)] until a self- 
consistent solution is reached. We find for a convergent 
solution, the q 2 order parameters are not completely in- 
dependent. In fact, these q 2 order parameters break up 
into q sets of q order parameters with the same magni- 
tude. Taking p/q = 1/3 or 1/4 as examples, their rela- 
tions are summarized in Table HI 
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TABLE I: Pairing order parameters for p/q = 1/3 (left) and 
p/q = l/4(right). a, b, c, and d denotes some complex num- 
bers depending on details, like the fermion density and U/t. 

These structures can be understood from the symme- 
try properties discussed above. The system is invariant 
under translation by one lattice site along x-direction and 
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(a): p/q=1/3 



(b): p/q=1/4 



FIG. 2: (Color online) Structure of the vortex lattice found 
from self-consistently solving the BCS mean-field Hamilto- 
nian, where ub = 1/3 for (a) and ub = 1/4 for (b). This is a 
contour plot of pairing order parameter A(r). The grey area 
means low superfluid density and locates the center of vor- 
tex cores. The intersection of vertical and horizontal straight 
lines indicates lattice sites. In both plots we use U = — 5.5t 
and rif = ni = 1/3 for (a) and = = 1/2 for (b). 



a simultaneous translation of k y by 2irp/q. Under this 



operation, A 



-» Al,, where j3' 



P + 1 (mod q) and 



V = 1 + 2 (mod q), thus these two order parameters must 
be equal up to a relative phase. To verify these relations, 
we show in Fig. [3](a) that our numerical solutions satisfy 
I iv = |A r rA z t|/(|A'||A r |) = 1 for V = I + 2(mod q), 
where T is a q x q matrix with Yij = ^+i( mo d q),j- This 
symmetry imposed relation works for any p/q, which im- 
plies that if A is non-zero, all A 2n ( mod q ^ are non-zero, 
i.e., zero and finite- momentum component must co-exist. 

The self-consistent solution also determines the relative 
phases. For p/q = 1/3, we find six degenerate solutions 
with (0i,02) = (±27r/3,±27r/3),(0,±27r/3),(±27r/3,0); 
for p/q = 1/4, we find = +tt/2 and either a, b ^ 0, 
c = <i = or c, d ^ 0, a = b = 0, therefore there are 
totally four degenerate solutions. One can see from Fig. 
[2] that this degeneracy can also be inferred naturally from 
the geometry of the vortex configuration. 

The most favorable relative phases can also be un- 
derstood by a simple Ginzburg-Landau (GL) argument. 
This GL theory should work well particularly nearby the 
phase transition point discussed below, where the order 
parameter is small. For those order parameters that defi- 
nitely co-exist, we first write down the most general cou- 
pling form between them by momentum conservation, 
and then determine the most favorable relative phases 
by minimizing energy. For instance, for p/q = 1/3, A , 
A 1 and A 2 all co-exist, and then one can write 

E GL oc A * A^A 1 A 2 + A 1 * A 1 * A A 2 + A 2 * A 2 * A^ 1 + c.c, 

thus the energy depends on the phases as cos(2#i — O2) + 
cos(2#2 — 0i) + cos(#i +62), and one can easily show that 
the angles listed above are its minima. For p/q = 1/4, 
A and A 2 definitely co-exist, thus one shall write down 

Eoc A°*A°*A 2 A 2 +c.c, 
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FIG. 3: (Color online) (a) For p/q = 1/3, I02 and hi (see 
text for definition) equal to unity within numerical accuracy, 
which verifies the symmetry relations, (b) Average superfluid 
order parameter A = Ai as a function of U for p/q = 1/3, 
n — 1/3 (red solid line) shows a phase transition and n = 1/2 
(blue dashed line) does not. 



which gives the energy-phase relative as cos2#i, whose 
minima occur at 6\ = ±7r/2. 

Vortex Configuration: To study the configuration of 
vortices in the superfluid ground state, we first note that 
in presence of magnetic field the Wannier wavefunction 
at each site should be chosen as 



i2irp(x/\)(y/\-j y )/q 



^(r-Rj), (9) 



where ip° is a Wannier wavefunction in absence of mag- 
netic field, A is the lattice spacing. The real space pro- 
file of the order parameter A(r) = ^ . Aj(/?(r — Rj) is 
contour-plotted in Fig. [2] for two different flux densi- 
ties. We also verify that the phase of A(r) winds 2tt 
around each vortex core. There are six (four) space group 
symmetry-related configurations for Fig. [2fa)(b), which 
corresponds to six (four) degenerate mean-field solutions. 
Hence, we have presented a systematical way to deter- 
mine the configuration of vortices in a BCS superfluid 
from a microscopic theory, which can be verified experi- 
mentally with standard imaging technique in cold atom 
experiments. 

Insulator (semi-metal) to superfluid transition: For 
= p/q, and for the fermion density of each spin com- 
ponent n = v/q, where v is an integer from 1, . . . , q — 1, 
the system is usually a Hofstadter insulator in absence of 
interactions. Except for the case that q is an even integer 
and n = 1/2, the system is a semi- metal since there are 
Dirac nodes at the Fermi energy. In both cases, since the 
Fermi energy is either in the band gap (Hofstadter insula- 
tor), or the density-of-state linearly vanishes (Hofstadter 
semi-metal) at the Fermi energy, there is no Cooper in- 
stability for infinitesimally small attractive interactions. 
Thus, it requires a critical interaction strength to turn the 
system into a paired superfluid through a second-order 
phase transition, as shown in Fig. [3f b). In this calcula- 
tion, we also fix the fermion density by judging chemical 



potential. This transition is driven by the competition 
between pairing-energy gain and the single-particle en- 
ergy cost to excite particles across the band gap, which 
was first discussed in Ref. [Icj] for a lattice system with- 
out magnetic field. Without the magnetic field, to realize 
the transition one needs to tune the interaction close to 
a Feshbach resonance to achieve strong pairing strength 
comparable to the band gap; while in this case, since the 
magnetic band gap is controlled by original band width 
t, the transition can be achieved by varying U/t, as rou- 
tinely done in cold-atom experiments. This transition is 
accomplished by a change in compressibility and can be 
measured directly from in situ density profile, which has 
been successfully used in studying boson Hubbard model. 
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